!-------------------------------------------------------------------------------------
! subroutine GetEfield
!-------------------------------------------------------------------------------------

subroutine GetEfield(dfu, kXvec, kZvec, ExFilter, phi02DFT, &   ! inputs 
                     Ez2D, Ex2D, OmegaBounce, Density2Dreal, &  ! outputs
                     density1DFT, phi2DFT, &                    ! outputs
                     CEz2D, EzKx)                               ! outputs


  use fft_wrap,   only : fft_forward, fft_backward
  use dimensions, only : jx, jz, jv, jDirection, dv, nkz
  use parameters, only : iunit

  implicit none



  !-- external variables --


  real,    intent(in)  ::  dfu(jz,jv,jx,jDirection) 
  real,    intent(in)  ::  kZvec(jz), kXvec(jx)
  real,    intent(in)  ::  ExFilter(jx)

  complex, intent(in)  ::  phi02DFT(jx,jz)

  real,    intent(out) ::  Ez2D(jx,jz)
  real,    intent(out) ::  Ex2D(jx,jz)

  real,    intent(out) ::  OmegaBounce
  real,    intent(out) ::  Density2Dreal(jx,jz)
  complex, intent(out) ::  density1DFT(jz)   ! diagnostics only (?)
  complex, intent(out) ::  phi2DFT(jx,jz)    ! diagnostics only
  complex, intent(out) ::  CEz2D(jx,jz)      ! diagnostics?
  complex, intent(out) ::  EzKx(jx,jz)       ! diagnostics?


  !-- internal variables --

  real     PhiCenter
  integer  ikx, ikz
  integer  ix, iz, iv, iDirection


  real     Phi1D(jx)
  complex  density1D(jz),  Phi1Dx(jx),  Phi1Dz(jz)
  complex  density2DFT(jx,jz)
  complex  density2D(jx,jz)

!-------------------------------------------------------------------------------------
      

  do ix=1,jx

     do iz=1,jz
        density1D(iz)=0.0
        do iDirection=1,2
           do iv=1,jv
              density1D(iz) = density1D(iz) + dfu(iz,iv,ix,iDirection)*dv 
           end do
        enddo
    end do !  end iz loop
    
    Density2Dreal(ix,:) = real( density1D(:) )
    
    call fft_forward(jz, density1D, density1D)
    density1D = density1D/float(jz)
    call dealias(jz,density1D)

    density1D(1)=0.0                  ! subtract uniform background
    density2D(ix,:) = density1D(:)    ! density2D(x,kz) 
                                      ! partial Fourier transform
                                      ! COMPLEX density2D(jx,jz), LOOKS OK
  enddo  !  end of ix loop

  do ikz=1,jz

     Phi1Dx(:) = density2D(:,ikz)       !    density2D(all x,kz)   complex Phi1Dx(jx)
     
     call fft_forward(jx, Phi1Dx, Phi1Dx)
     Phi1Dx=Phi1Dx/float(jx)
     call deAlias(jx, Phi1Dx)

     density2DFT(:,ikz) = Phi1Dx(:) * ExFilter(:)  !  density2DFT(all kx,kz), filtered in kx

     phi2DFT(:,ikz) = density2DFT(:,ikz) / (kXvec(:)**2 + kZvec(ikz)**2 )

  enddo

  phi2DFT = phi2DFT + phi02DFT          !!!!!!!!!    TOTAL POTENTIAL

  phi2DFT(1,:) = 0                      !  spatially uniform in x not allowed (03/21/2011)
        
  density1DFT(:) = density2D(jx/2,:)                !   diagnostic



  do ikx=1,jx                                       ! get Ez2D(x,z) complex CEz2D(jx,jz)
     Phi1Dz(:) = -iunit * kZvec(:) * phi2DFT(ikx,:)
     call fft_backward(jz, Phi1Dz, Phi1Dz)
     CEz2D(ikx,:) = Phi1Dz(:)                       ! Ez(kx, all z)
  end do

  EzKx = CEz2D                                       ! Ez(kx, all z),  needs  EzKx
  do iz=1,jz
     Phi1Dx(:) = CEz2D(:,iz)
     call fft_backward(jx, Phi1Dx, Phi1Dx)
     Ez2D(:,iz) = Phi1Dx(:)
  enddo	                                             !  got Ez(all x, all z)



  do ikx=1,jx   !  get  Ex2D(x,z)      
     Phi1Dz(:) = -iunit * kXvec(ikx) * phi2DFT(ikx,:)
     call fft_backward(jz, Phi1Dz, Phi1Dz)
     CEz2D(ikx,:) = Phi1Dz(:)                         ! Ex(kx, all z)  complex CEz2D(jx,jz)
  end do

  do iz=1,jz
     Phi1Dx(:) = CEz2D(:,iz)
     call fft_backward(jx, Phi1Dx, Phi1Dx)
     Ex2D(:,iz) = Phi1Dx(:)                   
  enddo                                               ! got Ex(all x, all z)

  density1D(:) = Ez2D(jx/2 + 1,:)

  call fft_forward(jz, density1D, density1D)
  density1D = density1D/float(jz)
  call dealias(jz, density1D)

  ikz=2

  ! amplitude of TOTAL potential at speckle center of sin(kz*z + phase)
  PhiCenter = 2 * abs(density1D(ikz) / kZvec(ikz))

  OmegaBounce = kZvec(ikz) * sqrt(PhiCenter)
     
  return

end subroutine GetEfield
    
    
